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Abstract 

We give a new fast method for evaluating sprectral approximations of non- 
linear polynomial functionals. We prove that the new algorithm is convergent if 
the functions considered are smooth enough, under a general assumption on the 
spectral eigenfunctions that turns out to be satisfied in many cases, including 
the Fourier and Hermite basis. 
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1 Introduction 

The goal of this paper is to introduce and analyze a new method to compute 

spectral approximations of polynomial functionals typically arising in spectral 
numerical methods applied to nonlinear partial differential equations. 

To describe the method and results, let us consider for instance the functional 

X{u){x) —uixf (1.1) 

where u(x) is a smooth function on the one-dimensional torus T and p > 2 an 
integer. We can expand u{x) as the Fourier series 

where the Ujt G C are the Fourier coefficients associated with u. In this case, the 
functional X{u){x) = J2k&^^'^^-^k{u) satisfies the convolution formula 

VfcGZ, Xk{u)= ^n---^jv (1-2) 

(ii,-jp)ezp 

k=ji-\ \-jp 



To compute a numerical approximation of such a quantity, a direct method would 
be prohibitive: if u is approximated by coefficients, the sum on the right-hand 
side involves N^^^ terms making the computational cost prohibitive for large N. 
That is why standard methods use the Fast Fourier Transform (FFT) algorithm 



to evaluate (1.1) on grid points and an inverse FFT to go back to approximated 
Fourier coefficients X^. Though this method has the disadvantage to introduce 
aliasing problems due to the structure of FFT, it is very cheap in the sense 
that if the grid is made of N points (and u approximated by N frequencies) the 
computational cost is of order log A'^. 

In many other situations like Hermite spectral methods, the problem is much 
harder because of the lack of fast transformation algorithm from collocation grid 
points to spectral variables (see however [10] for recent results by A. Iserles on a 
fast algorithm to compute Legendre coefficients). 

In this paper, we would like to show how a direct sparse approximation of 



( 1.2 ) of the form 



VfcGZ, Xf(n)= ^h---^Jr>- (1-3) 

k=ji-\ hip 

lii|-bpl<A^ 

yields a consistent approximation of in the sense that we can control the dif- 
ference — X^ll in some Banach algebra, provided the function u is sufficiently 
smooth. 



The big advantage of the representation (1.3) is that the sum on the right- 
hand side involves only C'(A^(log A^)^"^) terms making the direct approximation 
at the spectral level possible and efficient. 

To have an idea of why this method is valid, let us calculate directly the 
difference 

Xk{u) - X^{u) = ^ ^ji • • • 

fc=iiH — hip 
liil-bpl>^ 



^jp- 



We can write 



\Xk{u)-X^{u)\ < ^ \hnnn\---\fp\Wjp\ 

fc=iiH — hip 
lii|-lipl>A^ 



and we immediatly obtain the bound 
||X(n)-X»||^, ■.= Y\Xk{u)-X^{u)\ <^AY. 



fcez fcez 



iJll^-l 



P 1 II iiP 



(1.4) 

Note that here we used £"^-based spaces (Wiener algebras) as they are the sim- 
plest to deal with polynomial nonlinearities in spectral representation. Even if 
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similar results can be obtained for standard £^-based spaces, we will state our 
result in these Banach spaces to avoid too many technical details. As high reg- 
ularity i'^ and £^ spaces are imbricated, this does not affect the validity of our 
approximation results. 

We see however that the previous proof does not extend straightforwardly to 
the case of Hermite functions. In this case, if u{x) is defined on the real line M 
and decomposes into 

u{x) = "^UkXkix) 
fceN 

where the Xj(^)) J ^ are normalized Hermite functions, then the Hermite 



coefficients of the product (1.1) are given by 



VfcEN, Xkiu)= ak-ji,...,jpUj^---Ujp, (1.5) 

(ii,--,ip)eNp 

where 

= / Xk{x)Xji{x) ■ ■ ■ Xjp{x)dx (1.6) 
Jm. 

are the integrals of products of Hermite functions. Note that in this situation, 

the coefficients are non zero even in the case where k ^ ji + ■ ■ ■ + jp. To obtain 



a convergence result similar to (1.4) we thus see that we need a non trivial 
control of these coefficients. To this aim, we take advantage of the recent work 
by B. Grebert, R. Imekraz and E. Paturel, see 0, in which bounds are given 
for the coefficients afc:ji,...jj, that allow to prove that X acts on high-regularity 
Sobolev spaces. Note that this Hermite case is of particular importance because 
of the lack of fast Hermite transform, while in practice Hermite spectral methods 
are quite natural and widely used in many applications fields like Bose-Einstein 
condensate simulations and Fokker-Planck equations. 

In Section 2, we give a very general result in an abstract setting by assuming 



explicit bounds on the coefficients afc;ji,...jp in (1.5). This result covers the case 
of Fourier and Hermite basis, spherical harmonics functions, and eigenfunctions 
of operators of the form —A + V in dimension one with Dirichlet or periodic 
boundary conditions. 

To cover different situations, we introduce general sparse sets of indices of 
the form • • • \jp\ < N where a = or 1. 

In the case where the momentum k — ji — ■ ■ ■ jp is bounded in the sum defin- 



ing Xfc (like in the Fourier case, see (1.3)), the set of non zero coefficients will 



be indeed of size 0{N(\og Ny ^) for a = 0. However in more general situa- 
tions like Hermite approximation, the set in k and {ji, ■ ■ ■ , jd) will be of size 
0{N'^ {log N)P'^) for a = (if |A;| < N) and C'(A^(log iV)^) for a = l. The effect 
of this parameter a is only a slight deterioration of the rate of convergence of the 
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approximation, but it reduces drastically the computational cost of the method 
for large N in the Hermite case. 

In Section 3, we show how an iterative implementation of the algorithm yields 

a convergent approximation of the product of p functions with a cost of order 
0{pN log N) instead of C'(A''(log A'')^'"^). We give an error estimate for this case 
as well. 

In Section 4 we detail the case of periodic exponential functions (the Fourier 

basis) and discuss the possible extensions to cigcnfunctions of operators of the 
form —A + V. In Section 5 we consider the Hermite case and show by numerical 
experiments that the error bounds are optimal. 

2 An abstract result 

We consider Z = Z*^ or N*^ for d > 1. For u = {uj)j^z € we set 

ll«ll£i = Ellj'in^^-|' (2.1) 

where = max(l, . . . , |j'^|) for j = {j^, . . . G Z. We also define the 
norm 

\M,.={j:\\^t\uA')K (2.2) 

jez 

and using the Cauchy- Schwartz inequality, we can easily prove that if s' — s > d/2, 
there exists a constant C such that for all u, we have 

\\u\\i2 < WuWfi < C\\u\\^2 ■ (2.3) 

S S gf 

For a given integer j» > 2, we aim at approximating a function X : {ilY ~^ 
defined by X('u^, . . . , u^) = {Xg{u^, . . . , uP))£^z where 

yiez, Xe{u\...,un= Yl "'^m-iXr--^?.' (2-4) 

ji,-jpezp 

with given coefficients o,i;ji—jp G C We use the following notation: for a multi- 
index j = (ji, . . . ,jp) and £ & Z, we define the momentum 

M{e,j)=£-ji jp. (2.5) 

We will also sometime use the notation a^-j to denote the coefficient ae-jj^...jp. 
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2.1 Sparse sets of frequencies 

We consider a subset fC <Z Z. We will typically consider the case where /C = 
-E, a bounded set of or a sparse set of indices of Z. We assume that K, is 
equipped with a function | • | measuring the size of multi-indices of the form 
j = (j^, . . . , j'') G Z. We assume that there exist positive constants cq, Cq and 
a such that 

VjGZ, collill < |j| <Co||iir. (2.6) 
We then set for u G (compare 



W\,,=Y.\i\''\u,l (2.7) 
ie2 



and using (2.6) we obtain 

c\\u\\^,<\u\^,<C\\u\\^, (2.8) 

for some constant c and C independent of u. As particular cases of application, 
we mainly have in mind the two following situations: 

(i) /C is a set of the form 

ICM = {j eZ\\j\<M}, with |j|:=||j||, (2.9) 
where M G N can be equal to +oo in which case ICm = Z. In this situation, 



we have Cq = cq = a = 1 in the inequality (2.6). Note that for a given M, 
we have (J^a/ < CM'^ for some constant C independent on M, where ]\F 
denotes the cardinal of the set F. 

(ii) /C is a sparse set of the form 

d 

}C*M = {j eZ\\j\<M}, where |j| := J] (1 + |j-|), (2.10) 

n=l 

for some given M G N. In this case, using the inequality of arithmetic 



and geometric means, (2.6) is valid with a = d. In this situation, we have 
d^jv/ < CM(logM)'^~^ for some constant C independent of M (see for 
instance [H [TTj). 

For a fixed a G {0, 1} and > 0, we define the following approximation 

X^'"(n) = (Xf °),eyc oiX{u): 

V£G/C, Af"(ni,...,nP)= ^ a,,,,...,^v]^ ■ ■ ■ u^^. (2.11) 

l^nii|-bpl<JV 

The next Lemma estimates the number of non zero terms involved in the 
definition of X^'"(ti) in the two cases (i) and (ii) described above. 
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Lemma 2.1 The cardinals of the sparse sets of indices can be estimated as fol- 
lows: Let a G {0, 1} and p > 1. There exists a constant C depending only on d 
and p such that, for all M and N > 1, we have 



(i) With K-M defined by (2.9), and \j\ = \\j\\ for all j G Z, then 

UiJi, ■■■jp^ ICl, I l^rijil • • • <N}< C(tt/CM)('-"^iV"(logiV)P-^+", 

with the convention {\ilC]\j)^ = 1 when ICm = Z, that is M = +oo. 
(a) With /C^j and \j\ the sparse norm defined by (2.10), then we have 

m,ji,---jp(^inir\mji\---\jp\<N} 

Proof. The proof of (ii) is classical (see for instance dKH]) using the fact that 
in this case, \j\ = HLiCI + li'^l) when j = (ji, . . . so that 

d dp 

mh\ ■ ■ ■ \jp\ = ( 11(1 + i^'i))" n 11(1 + ij^i). 

fc=l k=ln=l 

which yields the result for a = 1 (independently on M). The case a = is 
treated similarly. 

The proof of (i) is a consequence of the fact that for all > 1 and p > 1, 

tt{ji, ■■■,jp€Zn lljill • • • llipll <N}< CpN''{logNy-\ 

for some constant Cp depending on p and d. We prove this by induction on 
p: for p = 1 the result is clear using ||j|| = max(l, |j^|, . . . , Ij"^!) G N\{0} for 
j = (j^, . . . G Z. Let us assume that it holds for p — 1 > 1. We have 

tJO'i,--- ,jp^2P\\\n\\---\\jp\\ <N} 

^ N 
= J2^{ji, ■ • • Jp-i G I lliill • • • IIVill < -} X Hi G Z I llill =k}, 

k=l 

N 
k=l 

^ 1 

< 2'^ dCp-iN'^ {log N)P-^ y^i.- CpN'^{\ogNY~^ 

k=l 

for some constant Cp depending on p and d. This yields the result. Here we used 
the fact that we calculate explicitly that for k >2, G 2^ | ||j|| = k} = dk'^~^, 
while for k = 1, this number is equal to 2"^, with the definition of \\j\\ . ■ 

As we will see below, the previous result can be refined when the coefficients 
aij in (2.4) have some special structure implying a decay property with respect 
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to the momentum defined in (2.5). We will consider theses cases more 

in detail in the section devoted to the Fourier case. 



2.2 Error estimate 

The goal of this section is to give an estimate of the error 
\\X{u\...,un-X'''''{u\...,unL, 



for smooth u^, . . . ,u''^ and where X^'" is defined by (2.11) for some given a G 
{0, 1} and > 1. 

We make a general hypothesis on the coefficients ai-j involved in the definition 
of the functional X. 

Definition 2.2 Let k = {ki, . . . ,kq) G Z'^ with q > 1 a multi-index. For n = 
l,...,q, we set iin{k) the n-th largest integer amongst \\ki\\ ,...,\\kq\\, so that 
we have fJ-i{k) > ^2{k) > fJ^sik) > • • • . 

We make the following hypothesis: 

Hypothesis 2.3 There exist z/ > 0, € [0, 1] such that for all R, there exists 
cr such that for all I ^ Z, and all j = {ji, . . . ,jp) E Z^, we have 



where k = = {iji, . . .,jp) 



(2.12) 



Let us make some comments on this definition. Such bounds (with 9 = 0) were 
used in several recent works O El [21 [H [3l [7] to prove long time existence results 
on nonlinear PDEs set on manifolds with different kind of boundary conditions 
(compact manifold, Dirichlet, etc.). It holds true in many situations where the 
a£j are products of the form ( |1.6| ) with functions Xk defining a Hilbert basis 
on a manifod M, like the Fourier basis on a torus. It is also valid (with 9 = 0) 
in the case of spherical harmonics, see [5l[6], and when Xk sue well localized with 
respect to the exponentials, see [Il[3] and Definition 5.3 of [7]. This last situation 
corresponds to the case where the Xk are eigenfunctions of an operator —A + V 
with Dirichlet boundary conditions in dimension 1, and with a smooth periodic 
potential V. 

More recently this was extended to Hermite functions basis diagonalizing the 
quantum harmonic oscillator operator, see [8]. In this case the previous bound 
holds true but for 9 = 1/2. 

The main result of this section is the following. 
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Theorem 2.4 Assume that the coefficients Ui-j of the function X(u^,...,uP) 
satisfy the Hypothesis 



2.3 for some constants u > and 9 G [0, 1], and let X 



be the approximation (2.11 ) defined for a € {0, 1}, > 1 and (/C, | • |) C {Z, 



satisfying (2.6) for some constant a > 1. Let k > d be fixed. Then for all s > 
and s' > max(cjs + 6k, (1 — 0)k + v), there exists a constant C such that for all 
N and for all functions tt* G i = 1, . . . ,p, we have the estimate 

\\X{u\ ...,uP)- X^'"(ni, ...,uP)\\,,< C7V-^(^'^') n |u%,^ , (2.13) 

i=i 

where 

Pis, s') = mill ( l^^l^ll^gi _ (1 _ _ ^\ . (2.14) 
V aa + 1 / 

To prove this Theorem, we will use the following technical Lemma. The proof 
of this Lemma is postponed to the Appendix. 



Lemma 2.5 Let k > d. Assume that ai-j satisfies the previous Hypothesis 2.3 
for some constants ly > and 9 S [0,1]. Then for all r > 0, there exists a 
constant Cr such that for all j = (ji, . . . ,jp) G Z^, 

E ll^ir M < am(jT+%2(j)(^-')'^+^ (2.15) 

Proof of Theorem \2.4\ We set for i e Z, 

Rt{u\...,un = X,{u\...,un-X^^''{u\...,uP) 



Knii|-bpl>JV 



For some t < s' , we can write 



Knii|-bpl>A^ 



Hence we get using (|2.6|) 

i=l 

where 

ih,--,jp)'^2P \Jl\ \Jp\ 
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where Cq is the constant appearing in (2.6). Applying the previous Lemma with 
r = as + aa{s' — t) and using again (2.6) we see that C{t) will be finite if 



or equivalently 



t = max(crs + aa{s' — t) + Ok, (1 — 6)k + i^). 



as + aas' + eHi a\ , \ 
t = max( , (1 — t))K + v) 



aa + 1 

in which case s' — t = /3(s, s'). This shows the result. 



3 Iterative approximations 

We consider now the case where X{u^, . . . , u^) corresponds to the product oper- 
ator of p functions = X^jg^ ''^)Xjix), 1 < ^ < where Xji^) is an orthonormal 
basis of L'^{M) where M is a manifold (typically M = T*^ or M'^). In this case, 
the coefficients a£;ji,...j„ are given by the integrals 



J M 



We will see in the example below that bound (2.12) holds in many situations 
such as the Fourier basis on T"^ and the Hermite basis on M^. In such a case, we 
identify a function u with its coefficients Uj and talk about n G £j by a slight 
abuse of notation. 

In the previous section, we have proven that for two functions and v?, 
the function yields a good approximation of the product u^u^ = 

X{v}-,v?) if these functions are smooth. Now for three functions v}, v? and n^, 
instead of approximating the product u^u'^u'^ by using u^, u^), which 

generates a computational cost of order 0(A^(log A^)'^) in dimension d = 1 and 



for a = 1 (see Lemma 2.1), we might use the following algorithm: 

1. Compute the approximation v = X^''^[v}' ,u^) of the product u^u^ 

2. Compute X^'"(t>, u^) as approximation of u^u^u^. 

In other words, we replace X^''^{v} ,v? ,v?) by A^'"(X^'"(n^, n^), u^). 

Obviously the cost of this algorithm is of order C'(2A^(log A^)'^) for a = 1, 
instead of 0(A^(log A^)^) (in dimension d = 1, see Lemma 2.1). Such an iterative 



approximation can be easily generalized to any product of p functions, and the 
global cost is of order C'(pA^(log Ar)2) for a = 1, instead of C'(Af(log A^)p). As 
we will see now, an error estimate of the same kind as in the previous section 
remains valid for such sparse approximations. For simplicity, we only present the 
result in the case where | • | = || • || , which implies \\u\\ ^ = \u\ ^. 
This is given by the following result: 
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Theorem 3.1 Let u^{x), i > be given functions. For all i > 0, let us define 
the functions J7^'"(n^, . . . ,u^) by induction as follows: U^'°^{u}) = and for 
i > 1, 

where a £ {0, 1} is fixed and X^'" defined in (2.11) for (/C, | • |) C {Z, || • || ) 
where \j\ = \\j\\ for all j G Z. Then for all p, the function U^'°'{u^, . . . ,vP) 
is an approximation of . . . u^) = ■ ■ - u^ in the following sense: Assume 

that the coefficients atj satisfy the Hypothesis 2.3 for some constants > 
and 9 £ [0, 1], and let k > d be fixed. Then for all p £ N, s > and s' > 
max(s + {p — 1)0k, (1 — 9)k + zv), there exists a constant C such that for all N 
we have the estimate 

\\X{u\ ...,un- U'''^{u\ ...,un\\,,< CN-P^^'^^'^ n 11^*11^,^ , (3.1) 

i=l 

where 

Pp{s, s') = min (^^^^-^^rr^' ^' ^^^'^ -i^-J^))- (3.2) 

Proof. As N and a are fixed, we set W := U^'°'{u^, . . . ,u^). For p = 2, the 
estimate is the one given in Theorem 2.4 with a = 1. Assume that it holds for 
p — l>2. In particular, we have for all s" > and s' > max(s" + {p — l)9n, (1 — 
e)K + u) 

p-i 

\\UP''\\,i <{l + CN-^'^-^^'"'''^)'[[\\u'\\^, , 

1=1 

for some constant C depending on s' , s" and p. Here we use the fact that in 
the case where \j\ = \\j\\ the norms || • || ^ and | • coincide. Now using the 
definition of U^, we can write 

UP - X{u\ ...,uP) = X^'"{UP-\uP) - UP-^ • uP 

+ {UP-^ - X{u^,...,uP-^))-uP. (3.3) 

As a direct consequence of Lemma [2. 5 [ we easily see that the following holds: for 
s > (1 — 26)k + u, and for u = J2jez "^jXj arid v = J2j£Z '^i^o ™ ^s; have 

\\uv\\.^< V \W \ai-j^j^\\uj^\\vj^\<Cs\\u\\^ 



pi— /, ll^ll \^^0\3l\\"-3\\\^32\ — ^sW^WfX ll'^IUi 



Using this inequality and (3.3) we obtain for s" > max(s + 6k, (1 — 6)k + i/) 
using (|3.1|) for p = 2, 



\UP -X{u\...,uP)\\^, < (:7iV-'^2(s,s")||[/P-i| 



\u 



p\ 



+ \\UP~'-X{u\...,uP-')\\^, \\uP\\^, 

s-\-Ok s-\-6k 
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and hence, for some constant C depending on s, s', s" and p, 



p 

i=l 

We take s" = s' - {p - 2)eK, so that /3p_i(s", s') = 0. For this s" we have 

P2\S,s ) = mm ,s + — — K 

V a + 1 

. (s' - s - {p-l)eK , 
= mm ^ , s — [p — 6)0k — k — v 



n I n l\ ■ ^ — S — Ul\, — \p — AUt\, I \ 

Pp-i[s + VK, s ) = mm I — , s — (p — A)&k — k — v]. 



Moreover, we have 

■ s' - s - Ok- {p-2)en , 
a + 1 

On taking the minimum between /3p_i(s + 6k, s') and /32(s,s"), we obtain the 
result. ■ 

In the rest of this paper, we will show how this Theorem can be applied to 
many situations including the discretization of polynomials in Fourier or Hermite 
basis. 

4 Fourier basis 

We consider now functions u{x) defined on x E T'^. We consider functionals of 
the form 

X{u^,--- ,vF){x) = h{x)u^{x)---vF{x), (4.1) 

where h{x) is a given function defined on the torus T^. With a function u{x) G C, 
x = (x^, • • • x'^) S T"^, and for a given j = [j^, ■ ■ ■ ^ Z := we associate 
the Fourier coefficients 



where j - x = j^x^ + • • • j'^x'^. In this case, the coefficients ai.j^...j^ defined in (2.4) 



can be calculated explicitely, and for given i ^ Z = 7L and j = (ji, . . . , jp) € Z'^ . 

^^^n-,. = E ^^7^^ /, e^(-^+^+^-^+-+^-^)-dx = (4.2) 
where the numbers hu are the Fourier coefficients associated with the function 



h{x), and with the definition (2.5) of the momentum A4{£,j). Here we use the 
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very special property of the exponential functions e^^'^ that the product of two 
basis functions is again a basis function. We assume that b{x) extends to an 
analytic function on a complex strip Up := T'^ x i[—p,p\'^ around the torus, 
which implies by standard Cauchy estimates that 

VfcGZ'^, < De-^l^l, (4.3) 

where D = sup^g^^ \b{z)\. 

With this calculation, we can prove the following result: 

Proposition 4.1 // the Fourier coefficients of the function b{x) satisfy the 



analytic estimate (4.3), then the coefficients a,ij^...j^ 
Hypothesis 2.3 with z/ = and = 0. 



defined in (4.2) satisfy the 



Hence we see that when b is analytic, we will have I3{s,s') 



^-XT in the 



formula (2.14), provided s' and s are large enough. The proof of the previous 



proposition can be found in [U [7] . As explained in these references, the same 
result holds true when the function u{x) is decomposed on a Hilbert basis ej{x), 
j G Z"^ that is well-localized with respect to the exponential. This includes in 
particular the case where ej(x) are the eigenfunctions of a differential operator of 
the form u{x) i— )■ —Au{x)+V{x)u{x) for some smooth periodic potential function 
V{x) in dimension d = 1. We refer to [7J for extensive discussions on the subject. 

Let us mention that in the particular case where b{x) is a trigonometric poly- 
nomial containing only a finite number of frequencies, the use of the parameter 
a = 1 is not mandatory to obtain sparse set of indices. This is a consequence of 
the Lemma below: 



Lemma 4.2 Considering the approximation (2.11), we assume that there exists 
q > such that 

\M{i,j)\>q^ae,j=0. 

The cardinals of the sparse sets of indices with a = can be estimated as follows: 
Let p > 1, then there exists a constant C depending only on d, q and p such that, 
for all M and N > 1, we have 



With Km defined by (2.9), and \j\ = for all j G Z, then 



■■■jp^ JCl, I IjiI • • • IjpI < iV} < CN''{logNr-\ 



(a) With IC*j^j and \j\ the sparse norm defined by ( 2.10[ ), then we have 

n^ju ■■■3v^ iniT I ■ ■ ■ \jp\ <m< civ(iogiv)*"^ 
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Figure 1: Sparse approximation of for a = 3. Left: convergence of the £^-error; 
Right: CPU time 



Proof. For a fixed 



in the sets considered, the estimates can be obtained 

Now under the 



similarly as in the proof of Lemma 2.1 , and are independent of • 
assumption on the momentum, if ai-j ^ then we have necessarily = 
m G Z'^ with \m\ < q. Summing in m then yields the result (with a constant 
proportional to q"^). ■ 

Note that the case considered in the introduction corresponds to b{x) = 1 
and g = in the previous Lemma. 

We show now on a numerical example the accuracy of the estimates above. 
We consider the function u(x) = Ylkez'^k^^'^^ with Uk = {1 + l^l)""^ so that 



u ^ i^, for s' < a — I. We compute by the direct method ( |2.11 ) and the 
iterative algorithm described in Section [3j In both cases we expect a maximal 



convergence rate 0(A^°+i ) in £ (that is for s = 0). In figure [ij (left) we plot in 
log scale the £^-error versus the sparse level in the case p = 3 for the different 
approximation methods. In figure [T] (right) we plot the estimated CPU time 
together with the theoretical bounds CiV(log A^)^"^"*"*^ for the direct method and 
CpN (log Ny^" for the iterative one. For convenience we plot only the theoretical 
bounds for a = 0. The version a = is clearly more accurate than a = 1. On 
the other hand it has only a minimal extra cost, so for this particular example 
it is clearly preferable. It should be pointed out, however, that this is due to the 
very simple form of the functional X{u) = for which ai-j = if \M.{£,j)\ ^ 0. 

Figure [2] shows the error versus the CPU time. It is clear from this plot the 
advantage of the iterative algorithm with respect to the direct method, as well 
as the advantage of a = with respect to a = 1. 

Finally, in figure [s] we show the convergence of the £^-error, still in the case 
p = 3 but for different values of a. We consider here only the case of a = 1 and 



the direct formula (2.11). The results in the other cases are analogous. For all 
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Figure 2: Sparse approximation of for a = 3. ^^-error versus CPU time 



computation of u , a=1 , direct approx. 




Figure 3: Convergence of the sparse approximation of with the direct formula (2.11 ) 
and a = 1 



values of a we recover the expected theoretical rate of convergence. 



5 Hermite 

We consider now the case where u{x) is defined on the real line (x G M) and the 
basis (xj)jeN is given by the set of normalized Hermite functions defined by the 
formula 

Txj ■■= -^(^) + ^^Xiix) = (2j + l)xA^), j G N, (5.1) 

with the condition ||Y1|Lo,^^ = 1- Note that here, with the notation of the 
previous sections, we have Z = N. For all j G N, the Hermite functions are given 
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by 



where Hn{x) is the n-th Hermite polynomial with respect with the weight e 
Recall that these Hermite polynomials satisfy 



EN, / Hn{x)Hm{x)e-''''dx = 2''n\y/^5r, 
Jr 



and the induction relations: 

Ho{x) = 1, and Hn+i{x) = 2xHn{x) — 2nHn-i{x), n> 1. 

In this situation, {xj{x))jeN is a Hilbert basis of L^(M) and for a given real 
function u{x), we can write 



"^UjXjix), where uj = u{x)xjix)dx. 
Here, note that the Hilbert space associated with the norm defined in 



u(x) = 



(2.2) coincides with the domain of the operator T** (see (5.1)). Using standard 
notations, the classical space defined by 

= {u{x) G H'{R) I X ^ x^d^uix) G L'^{R) for 0<p + q<s} 

corresponds with the domain of the operator T'^/^ (see for instance [9j) and hence 
with particular, we can write owing to (2.3) 

c||«|| < Ikll^i^^ < C\\u\\ g^, , 

provided s — s' > d, and for some positive constants c and C independent of u. 

Let us now consider the functional X{u){x) = u{xY for p G N. In this case, 
the coefficients a£^ji-jp in (2.4) are given by the formula, for (£, ji, . . . ,jp) G N^"'"^. 

= / Xe{x)Xji (x) ■ ■ ■ Xjp{x)dx. (5.2) 
Jr 

The following Proposition can be found in [8j, Proposition 3.6: 

Lemma 5.1 For all v > 1/8 and all R > 0, there exists cr such that for all 
£ EN, and all j = (ji, . . . ,jp) G W, we have 

IX^jkY ( /.2(fc)i/X3(fc)i 

\ai-3\<CR ^( Y Y ) . (5.3) 

/il(fc)24 ^/i2(fc)2/i3(fc)2 +/il(fc) - /i2(A;)^ 

where k = = {£,ji, ■ ■ ■ ,jp)- In particular, these coefficients satisfy Hypoth- 

esis 



2.3 with e = 1/2 and v > 1/8. 
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Note that the estimate given in [8] is shghtly better than the one given in 



2.12 



because of the presence of the term /ii(fc)24 in the denominator in (5.3). 



Remark 5.2 In this paper, we will only consider the case of Hermite functions 
in dimension 1. The extension to higher dimension can he made using the frame- 
work of I3i, Section 3.2. 



N- 



In this situation, and when a 

\ — s'-\-K 



2.4 



1, we obtain a convergence rate of order 
2 for K > 1 and sufficiently large s' for the algorithm descr ibed in Theorem 



and 



-s' + {p-1)k/2 



for the iterative algorithm (see Theorem 
In the following we will illustrate these results by numerical 



3.1) 



simulation. In 

all the computations presented below, the coefficients \b.2\ are approximated 



in double machine precision by using Gauss-Hermite quadrature rules with the 
packages provided by J. Burkhardlj^ 

We first consider the case where p = 3, and for given number a, we consider 
the fonctions u{x) = J2n>o'^riXn{x) with m„ = (1 + n)""" so that u G il, for 
s' < a — 1. Hence in this case, we expect a maximal convergence rate of order 




Figure 4: Convergence of the sparse approximation 

0{N 2 ) in £^ (that is for s = 0) and when a = 1. In Figure[4]we plot in log- 
log scale the error measured in norm between the approximation X^'^{u, u, u) 

^ http: / / people.sc.fsu.edu/ ~jburkardt / cpp_src /hermite _rule /hermite jrule.html 



16 



and the exact solution v? whose Hermite coefficients are approximated using a 
Hermite transform with 500 points. The convergence rates observed correspond 
to the theoretical estimate (2.13). 

In figure [5} we plot the time required by the algorithm in the cases p = 2, 3 
and p = 4, to compute the Hermite coefficients of u{x)p. As expected, the time 
increases when p becomes large, which is in accordance with the cost of order 
A^(log A^)^ predicted by Lemma 2.1 We compare with the cost of the Hermite 
transform algorithm with N points, which is in ©(iV^). Note that for this latter 
method, the cost does not significantly differ with p, and only p = 3 is shown. 



0.014 




Sparse level N 



Figure 5: Time versus sparse level N 

In figure [6| we give the same time computation but with the iterative algo- 
rithm. We observe a significant speed up in the algorithm in comparison with 
the previous algorithm. 

In the last figures [7| [s] and [oj we fix a = 10 for the coefficients Un = (1 + n)~'^ 
of the function u{x), and we plot the error versus the time required for the 
algorithm (obtained in Figure [s]). We compare with the result obtained with 
the Hermite transform method. The results obtained are better for the sparse 
approximation. The results obtained for the iterative algorithm are similar, but 
less convincing because it requires much large number N, despite a better cost 
for a single iteration. 



17 




Sparse level N 



Figure 6: Time versus sparse level (iterative algorithm) 



Appendix: Proof of Lemma 2.5 



To prove this Lemma, we will use the following result, which can be found in ^ 
for 6* = and [8j for 6 = 1/2. 



Lemma 5.3 Assume that j = (ji, . . . ,jp) G Z'^ and for i £ Z, let 



Then we have 



M2(^,j)V3(^,jr-'+/"l(^,j)-/^2(Aj) 

yiGZ, \\e\\Ae{i,j)<2f,i{j). 



(5.4) 



(5.5) 



Proof. If ll^ll < the equation (5.5) is obvious using the relation Ag{£,j) < 

1. 

In the case \\£\\ > /ii(jf), then we have /^2(^, j)^/"3(^, j)^~^ = /"i(j)^Ai2(j)^~^ and 



m(j)V2(j) 



1-6 



W(j)V2(j)l-^ + ril -/il(j) 
ll^ll-m(j) \ 

M3r^^2{jy-' + ¥\\-Mj) 



m(j)V2(jr-'+m(j)^e(^,j). 




Figure 7: a = 10, p = 2 



We conclude using tlie fact that A0{£,j) < 1, and 1 < /Ui(j)^/i2(j)^ ^ < 



2.5 



Proof of Lemma 

three cases in the sum in ^ E ^ appearing in (|2.15l). 



Let j = {ji, ■ . . ,jp) G be fixed. We distinguish 



(i) > fJ,i(j). In this case, we have = \\i\\, fi2{i,j) = ^i(j) and 

l^si^jj) = /^2(j)- Hence we can write using ( 2.12[ ) with R = r + k, and the 
previous Lemma 

\W\>M3) Pll>m(i) 

< c2>i(j)>2(jr E M^jr- 

m\>Mj) 



As in this case we have 

Aii(j)V2(j)'- 



< m(j)V2(j)'" 



we obtain 

IKII>Ml(j) IKII>Mi(j) " " 



where the constant C depends only on r and k but not on j. Here we use the 
fact that the sum in the right-hand side is convergent and independent of 
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° Hermite transform 
— Sparse approximation 




x10 



Figure 8: a = 10, p = 3 



owing to tiie condition k> d. 

(ii) ^i(j) > > /i2(j). In this case we have = /U2(^, j) 

and lJ-3{i,j) = ^2(j)- Hence we have 



1-6 



1 



1 + 

Using again the previous Lemma, we obtain 

MlO)>PII>M2(i) 

1 



E 

w(i)>IKII>/^2{i) 



MiW)>IKII>M2W) 

and we conclude as in the previous case. 

(iii) /i2(j) > ll^ll • In this last situation, we can estimate directly the term, and 
obtain using the fact that A{£,j) < 1, 

E ii^iri«^;ii</^2(jr E i«^;ii<c/i2(jr+'^( E i)- 

IKII<M2(i) ll^ll<M2(i) IKII<M2(i) 
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Figure 9: a = 10, p = 4 
Hence, using ^{l E \ < ^2(j)} < C/X2(j)^ we get 

m\<f^2U) 

as K > d. Gathering the previous estimate yields the result. 

■ 
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